Data for testing IHS

library(tidyverse)
library(terra)
library(sf)
library(leaflet)

Obtain DEM

Pulling out a patch of publicly available LiDAR near Levin, Manawatū-Whanganui, on New Zealand’s North Island, from the LINZ data portal.

This dataset can help me identify the tiles I want. This has been downloaded as a geopackage and unzipped to data_spatial/aoi.

tiles <- 
  read_sf(file.path('data_spatial', 'aoi',
                    'manawatu-whanganui-lidar-1m-index-tiles-2015-2016.gpkg'))

# had a gander in QGIS, decided I want this bit:
wishlist <- tiles %>% 
  dplyr::filter(TileName %in% 
                  c('BN33_1000_1444', 'BN33_1000_1445',
                    'BN33_1000_1544', 'BN33_1000_1545'))

aoi <- st_as_sf(st_union(wishlist))

The data itself is downloadable in GeoTiff format from this dataset.

Unfortunately the dataset only has view services (WMTS and XYZ) and a point-data query API, so the tiles have to be downloaded manually from the ‘tiles table’ interface, e.g.:

This has been done and the four tiles unzipped and placed in data_spatial\elevation

elev_dir <- file.path('data_spatial', 'elevation')
list.files(elev_dir, pattern = 'DEM.*\\.tif$')
#> [1] "DEM_BN33_2015_1000_1444.tif" "DEM_BN33_2015_1000_1445.tif"
#> [3] "DEM_BN33_2015_1000_1544.tif" "DEM_BN33_2015_1000_1545.tif"

The easiest way to start working with the tiles is to create a VRT over them using GDAL. Below, my OSGeo4W GDAL installation is accessed to do this.

if(!file.exists(file.path('helpers', 'tile_files.txt'))) {
  tile_files <- 
    list.files(elev_dir, pattern = 'DEM.*\\.tif$', full.names = TRUE) 
  txt_con <- file(file.path('helpers', 'tile_files.txt'))
  writeLines(tile_files, txt_con)
  close(txt_con)
}

source(file.path('helpers', 'gdal_env.R')) # environment variables

system2('gdalbuildvrt',
        args =
          c('-input_file_list', file.path('helpers', 'tile_files.txt'),
            file.path(getwd(), elev_dir, 'lidar_1m.vrt')))

This can now be opened in R just like any other raster data source:

dem_file <- file.path(getwd(), elev_dir, 'lidar_1m.vrt')
dem <- terra::rast(dem_file)
dem
#> class       : SpatRaster 
#> dimensions  : 1440, 960, 1  (nrow, ncol, nlyr)
#> resolution  : 1, 1  (x, y)
#> extent      : 1800640, 1801600, 5503200, 5504640  (xmin, xmax, ymin, ymax)
#> coord. ref. : NZGD2000 / New Zealand Transverse Mercator 2000 (EPSG:2193) 
#> source      : lidar_1m.vrt 
#> name        : lidar_1m

A hillshade would also come in handy,

system2('gdaldem',
        args =
          c('hillshade', '-multidirectional', 
            dem_file, 
            file.path(getwd(), elev_dir, 'hillshade_1m.tif'))
        )

hs <- rast(file.path(getwd(), elev_dir, 'hillshade_1m.tif'))

---